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AN APPROXIMATE METHOD FOR SOLVING 


THE SINGLE -DEGREE -OF -FREEDOM ROLL EQUATION 

WITH TIME -DEPENDENT COEFFICIENTS 

By C. William Martz 
Langley Research Center 


SUMMARY 

/S033 — 

An approximate analytical method (with three alternate integration pro- 
cedures) is developed for solving the single-degree-of-freedom roll equation 
with time-dependent coefficients. The method is applied with each integration 
procedure to two sample problems and shows good agreement with more exact 
numerical solutions. Information governing the approximate error of each inte- 
gration procedure is also presented. 


INTRODUCTION 



One method often used to improve the accuracy of predicting the flight 
trajectories of rocket vehicle systems has been that of spinning the vehicle 
slowly very soon after launching. This tends to average out the effects of 
any asymmetrical vehicle loads on the trajectory. In other instances a vehicle 
may be designed to spin at higher rates for the purpose of spin stabilization. 
In either case the spin is usually provided in flight by means of deflected 
control surfaces, canted or twisted fins, or with small spinner rockets or gas 
jets . 


During the preliminary design of these vehicles, the fin cant (if this is 
the method chosen to produce the rolling moments) must be determined for each 
stage to produce the desired spin rates. This computation involves solutions 
of a first-order differential equation having time-variable coefficients. 

These solutions are carried out with a step-by-step numerical integration pro- 
cedure usually with the aid of a high-speed digital computer. 

The present paper presents an approximate analytical method (with three 
alternate integration procedures) for solving the single-degree-of-freedom roll 
equation with typically varying coefficients. The use of this method avoids 
the tedious step-by-step method of solution and allows quick and accurate hand 
computation of the results. The method is applied to two actual flight vehicle 


situations and the results compared with more exact solution determined 
numerically.^ 

In this paper, the approximate method is applied only to the roll equa- 
tion. However, the method can be (and has been) applied successfully to other 
physical problems. 


SYMBOLS 

The units used for the physical quantities defined in this paper are given 
in both the U.S. Customary Units and the International System of Units (Si). 
Factors relating the two systems are given in reference 1. 

A,a,D 

constants for simulating roll damping terms, second - ! ( see 
eq. (4)) 

Aj_ = Ae a t + 

D 

Aq = Ae & ! + 

D 

B,b,C,c 

constants for simulating roll driving terms, second"! (see 
eq. (4)) 

c l, c 2 

proportionality factors, seconds per radian and radians, 
respectively 

Ei,E 2 

complex Fresnel integrals defined by equations (9 a ) 


inertia of vehicle about spin axis, slug-foot^ (kilogram-meter^) 

ki,k2,k 5 

constants defined by equations (l5a) to (l5 c ) 

M X,P 

roll damping moment per unit p, foot-pound-seconds per radian 
( joule -seconds per radian) 

%,6 

aerodynamic roll driving moment per unit 6, foot-pounds per 
degree (joules per degree) 

M j 

jet damping moment coefficient, foot-pound-seconds per radian 
(joule-seconds per radian) 

M x 

time-dependent roll input moment, foot-pounds (joules) 

P 

vehicle spin rate, radians per second 


■^The more exact numerical method uses a fourth-order Runge-Kutta procedure 
with extrapolation to zero interval size as a correction factor and with tabu- 
lar look-ups for the coefficients. 
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Ap approximate error in spin rate, radians per second 

p spin rate at midinterval time, radians per second 

P 0 spin rate at zero interval time, radians per second (in multi- 

step problems, this quantity is final spin rate of previous 
interval. ) 

p ss steady-state spin rate, radians per second 

t interval time, seconds (t = 0 at beginning of each interval) 

t midinterval time, seconds 

u independent variable used in complex Fresnel integral 

y transformation variable defined by equation (5) 

S differential deflection of control surfaces or fins producing 

rolling moment, degrees 

Ae a t error in e a "^ (see eq. (13)) 


Subscript : 

o zero interval time 

A dot over a symbol indicates the first derivative of the quantity with 
respect to time. 


ANALYSIS 


The single degree-of-freedom rolling-moment equation can be written as 
follows : 


!XP + 

%,pP 

= MXjb 5 

+ M x 

+ Mjp 

(1) 

Inertial 

acceleration 

moment 

Aerodynamic 

damping 

moment 

Aerodynami c 

driving 

moment 

Auxiliary 

driving 

moment 

(time 

dependent ) 

Jet 

damping 

moment 



The aerodynamic driving moment can be generated by deflected control surfaces, 
canted fins , or any asymmetric rolling -moment condition dependent upon aero- 
dynamics. The M x term was included for time -dependent rolling moments such 
as would be obtained with small spinner rockets or gas jets. The jet damping 
moment, included for completeness, is rather hard to predetermine and is highly 
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dependent upon the configuration of the rocket motor propellant. In the fol- 
lowing analysis, the jet damping term is considered absorbed within the aerody- ' 
namic damping term. For conditions where the coefficients in equation (l) are 
constants, the solution can be written as 

P = (p 0 - Pss) e Ix + Pss ( 2) 


where p Q is the initial spin rate and the steady-state spin rate is 

Mx, 5 5 + M x 


Of primary interest in this paper is the condition in which the coeffi- 
cients of equation (l), normalized with respect to Ix, are not constant. The 
normalized equation to be solved is 


P + 


ix 


^.S 5 + Mx 

ix ix 


(3) 


In this equation, Ix can change continuously during a thrust period. The 
terms %,p and Mx^gS are both functions of air density, vehicle velocity, 
and Mach number which change throughout a flight; also, M x /lx must be given 

the capability of varying with time. The general form of equation ( 3 ) , chosen 
from practical considerations to simulate these possible variations, is 


p + p^D + Ae a "k) = Be^ + Ce 0 ^ (4) 

The method for solving equation ( 3 ), then, consists of fitting actual problem 
histories of the normalized driving and damping moments to the constant and/or 
exponential coefficients of equation (4) which can be integrated for spin rate. 

It is not implied that g&^Ix I s simulated by Be^ and M^I^ by 

CeCt. Rather, it is suggested for reasons of flexibility that the sum of the 
input terms be simulated by the sum of the two exponentials . Since the two 
input terms on the right of equation (4) are of identical form, their partic- 
ular solutions will be of the same form. Thus, for purposes of solving equa- 
tion (4), only the B term will be considered. The response of the C term 
will be added to this solution by the principle of superposition. 

Three procedures are presented for integrating equation (4). These pro- 
cedures referred to as the asymptotic method, the tabular method, and the 
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mid-damping method are discussed in the following sections. With all three 
methods, it is sometimes necessary to divide total problem time into two or 
more steps or time intervals as explained in the section entitled "Limitations 
and Applications." The method to be used is then applied once in each time 
interval, the final value of spin rate in an interval being used as the initial 
value of spin rate for the next interval. 


Asymptotic Method 

In the asymptotic method, the dependent variable is first transformed as 
follows. If 


^e at +Dt 

y = pe a 

equation (4) (neglecting the C term) becomes 


(5) 


bt+^e a ^‘+Dt 
y = Be a 


( 6 ) 


This equation can be integrated repeatedly by parts in two different ways. 
Each result is in the form of an asymptotic series which transforms through 
equation (5) to a separate spin rate solution. The more useful solution is 
the following convergent series: 


Be bt 


b + D 


A 

1 ££ + 


(Aeat) £ 


(A e at h 


a + B + D (a + b + D)(2a + B + D) (a + b + D)(2a + b + D)(3a + b + D) 

A. . A 2 


- |(eat_i)_Dt J B 

+ e a <P, 


b + D 


1 - 


+ b + D (a + b + D)(2a + b + D) 


A^ 


(a + b + D) (2a + b + D)(3a + b + D) 


+ similar C terms 


(7a) 


In application, each series should he extended until the first term not used 
is "negligible." 

The second solution, an asymptotic series useful for very small values of 
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Be^t 


Ae‘ 


at 


1 - b + D - a + (b + D - a)(b + D - 2a) _ (b + D - a)(b + D - 2a) (b + D - 3a) + 


Ae a b 


(Aeat) £ 


( Ae at)3 


+ e 


- |(e a t-l)-Dt J b 
\ P ° ' A 


b + D - a + (b + D - a)(b + D - 2a) 
A A 2 


(b + D - a) (b + D - 2a)(b + D - 3a) 


+ similar C terms 


(7b) 


As before, the series are extended -until the first terms not used are 
"negligible 

Normally, equation (ja.) is a fast, efficient, and accurate means of 
approximating the solution of equation ( 3 ). However, if convergence is too 
slow, another of the integration methods can be used. 


Tabular Method 

The tabular method picks up the integration of equation (4) at equa- 
tion (6) which is 

bt+§e at +Dt 
y = Be a 


Approximating e a t by 1 + at + rearranging terms, and integrating give 

2 

A A _ (A+b+D) 2 rt Aa / , A+b+p f 

y = p e a + Be a 2Aa / e 2 \ Aa ' dt 
0 2 0 


The independent variable is transformed by letting 


2 

Aa A ,A + b + D\ _ . it 

TV ~ka ) = 1 2 


u c 


Then 
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The integral is now in the form of the complex Fresnel integral tabulated in 
reference 2. 


The p solution with the C term included by superposition is 


P 


(A+b+D) 2 
2Aa ifirt 





Complex Fresnel integrals 
tabulated in reference 2 


( 8 ) 


(9a) 



The spin acceleration p can be determined from equation (4). Note that this 
solution does not degenerate for zero damping (A = D = 0) or constant damping 
(a = 0). However, for A = D = 0 a straightforward integration of equation (4) 
yields the following solution: 
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and, for a = D = 0, 


B 


C 


P 



+ 


A + b 




+ 


A + c 



p = -Ap + Be^ + Ce c ^ 



Constant 

damping 

coefficient 


J 


( 10 ) 


Mid-Damping Method 

The use of a constant damping coefficient in equations (10) introduces the 
third method for integrating equation (4) - that is, the mid-damping method. 

In this method the value of the damping coefficient at the midinterval time is 
used throughout the interval and the B and C terms are exponential as 
before. Therefore, the terms D + Ae a t are replaced with their midinterval 
value Ap and the constant damping solution becomes 

p = p 0 e“^-l^ + — — ^ (e^ - e + — — - Let _ e j (ll) 

Aq + b V ’ Aq + c \ ' 

where 


Aq = Ae at + D 


LIMITATIONS AND APPLICATION 


The method for solving the roll equation consists of fitting actual his- 
tories of the driving and damping moments, normalized with respect to Ix> to 
the constant and/or exponential coefficients of equation (4). The equation 
is then integrated in one of three ways. This curve fitting is a limitation 
on the accuracy of the method. However, the errors generated depend upon the 
skill of the curve fitter and in any case can be kept as small as desired with 
the use of more intervals. When fitting the coefficients of the damping term, 
first plot the values of Mx^p^Ix on semilog paper as a function of time. If 

the plot is essentially a straight line over a given time interval, the constant 
term D is zero and the values of A and a are easily determined from the 
zero intercept and slope of the faired curve. If the plot is curved, determine 
by trial and error the value of D required to straighten the plot over the 
required time interval and compute the values of A and a from the straight- 
ened and faired curve. 
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Think now about fitting the normalized input moments of equation (3) by- 
considering a rolling-moment input which increases linearly with time from 
2 ft-lb (2.72 J) to 6 ft -lb (8.l6 J) over a 10-sec period. One exponential 
term cannot adequately simulate this curve. However, by thinking of fitting 
an exponential curve having a displaced origin, one can visualize a flatter 
curve which will more closely simulate the requirements as shown in figure 1. 



Figure 1.- Exponential simulation of a linear moment variation. 


This, of course, is the same as using two exponentials - one with a zero 
exponent - to meet the need. However, this improved simulation requires cal- 
culations of increased quality and quantity since the additional term involves 
the difference of two large numbers . 

In addition to curve-fitting errors, the accuracy of the method is limited 
by approximations which were used in the integration procedures. Consider first 
the tabular method in which the approximation 


e at ~ 1 + at + 

2 
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was used. As shown in the appendix, this approximation causes errors in spin 
rate over an interval as given approximately by 


Ap = Be*)t + Ce ct - p(Ae a t + p) ^ at 

ae a ^ 

where 

Ae at = e at - ^1 + at + 

The absolute value of the spin-rate error becomes larger as the value of at 
is increased with the result that solutions usually must be divided into two 
or more steps or time intervals. Equation (12) can be used to determine these 
step lengths for a given accuracy requirement. However, a much faster and more 
desirable method is to use step lengths consistent with a factor-of -2 change in 
the damping term This ” factor-of-2 criterion” is consistent with good 

accuracy (about 2 percent error). If more or less accuracy is desired, a 
smaller or larger factor can be used and the relative error will be governed by 
equation ( 13 ). 

Now consider the mid-damping method. The use of a constant (midinterval) 
damping coefficient over the entire interval causes a loss in accuracy toward 
the midpoint of the interval. The approximate value of the spin-rate error, 
developed in the appendix, is given by 


( 12 ) 


(13) 


Ap « k]^t + — 


k 2 


- a 


(a-Ai)t 
e ' x/ - 1 


^_^(a+b)t - lj - D J (p - p) dt (14) 


where 


tl - k 2 e( a ' Al 4 + k 3 e< a+b > t 


(15a) 


k 2 = A Po " 


_B 

Ai + b 


(15b) 
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k 3 = - 


AB 


Aj_ + b 


(15c) 


fo ( P ' ’) 


B, 


dt = -pt + — 


; / A l 


(Ai + b 



e~ A 2_t _ A 


B/b 

Aq + b 


(• 


,bt 


(I5d) 


Equation ( it) shows that the absolute value of the error builds up near the 
center of the intervals while remaining lower near the end points. Also, the 
use of longer intervals does not necessarily decrease the accuracy of the 
method near the end points. Thus, longer intervals can be used accurately if 
only the final spin rate of each interval is required. 

Finally, with the asymptotic method of integration, no additional errors 
are involved since approximations are not used with this method. 

In some problems, the semi log plots of My^p/lx, My f ^8 , or M x jl^ have 

discontinuities where the rocket motor cases drop off or torque motors start or 
stop, and so forth. These discontinuities are the usual cause for additional 
time intervals, and are illustrated in figures 2 and 3> the data plots for the 
sample problem. 


Sample Problem 

The approximate method was used separately with each integration procedure 
to simulate the flight spin history of the second stage of the Trailblazer II 
vehicle described in reference 3* These simulations start with the separation 
of a first-stage rocket from the second-stage configuration during an exiting 
trajectory. Ignition of the second-stage rocket motor occurs at this separa- 
tion and the motor thrusts about 6 sec followed by a 20-sec coasting period. 
During this 26-sec post separation period, the second-stage configuration is 
exiting the sensible atmosphere, and its spin rate must be increased from about 
0 to 65 rad/sec by means of precanted booster fins. The problem is to prede- 
termine the fin cant required to produce this increase in spin rate and in the 
process to generate the spin history of the vehicle over the 26 sec. 

It is assumed that a particle trajectory has been computed to furnish time 
histories of Mach number, velocity, and dynamic pressure for converting the 
aerodynamic driving and damping coefficients into the time histories of M^p^I; 

and , 5 /lx shown in the semilog plots of figures 2 and 3 7 respectively. 

These plots were faired with straight-line segments for reasons previously 
discussed. For these straight-line variations, the C and D terms of equa- 
tion (4) are not needed and are set equal to zero. Both plots have natural 
breaks at the end of thrusting (41.2 sec), and the break at 48 sec was chosen 
to best fit the curves . 
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Figure 2.- History of the roll damping term for Trailblazer II second-stage flight. 
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Flight time, sec 


Figure 3.- Variation of 


%,5 

ix 


during Trailblazer II second-stage flight. 
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In application of the tabular method to the problem, intervals were 
determined by the factor-of-2 criterion and modified (at no loss in accuracy) 
with equation (12) to the intervals given in table I. 


CONSTANTS FOR PROBLEMS 1 AND 2 


Problem 1 : 


Constants 

Step 

1 

2 

3 

4 

5 

A, sec -1 ...... 

0.244 

0.114 

0.0433 

0.0146 

0.00296 

b/ 5 , sec"l/deg . . . 

8.85 

6.65 

2.48 

0.82 

0.162 

C, sec“l 

0 

0 

0 

0 

0 

D, sec -1 

0 

0 

0 

0 

0 

a, sec“l 

-0.1359 

-0.3022 

-0.3022 

-0.2280 

-0.2280 

b, sec -1 

-0.0621 

-0.3078 

-0 . 3078 

-0.2317 

-0.2317 

c, sec"l 

0 

0 

0 

0 

0 

t, sec . . 

2.8 

1.6 

1.8 

3.5 

3.5 

t 0 , sec 

0 

5.6 

8.8 

12.4 

19.4 

Initial p Q = 0 rad/sec 


Problem 2: 

Constants for problem 2 are the same as for problem 1 except initial 

P 0 = 200 rad./ sec . 

Although these time intervals are tailored to the tabular method, they were 
used also with the mid-damping method for purposes of comparison. Equa- 
tions (8) and (ll), with constants determined from figures 2 and 3 and given in 
table I, were used to generate the spin rate histories shown in figure 1+ for 
these two approximation methods and compared with the more accurate results of 
the numerical integration method. Continuous curves are illustrated for these 
two approximate methods. In practice, however, only end points of the inter- 
vals would be computed with the mid-damping method. These end points would 
then be faired for the final spin-rate history. Also shown in figure 4 are 
several values computed by the asymptotic method (eq. (7))- The use of the 
asymptotic method reduced the number of problem intervals to three - each cor- 
responding to one of the straight-line segments of the semilog plots in fig- 
ures 2 and 3 and to steps 1, 2, and 4 of table I. All results were computed 
for 5 = 1.75° , a value that was determined by trial and error to produce the 
required spin rate. Since spin rate at any given time is proportional to S 

l 4 
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O 

d) 

(/) 

TJ 

03 

M 


03 

03 

M 

C 

•H 

a 
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60 

50 

40 

30 

20 




Method 

* Numerical 

Mid-damping (eq*(ll)) 

Tabular (eq.(8)) 

O Asymptotic (eq.(7)) 



J I i I I — 

8 12 16 20 24 


28 


Time from separation, sec 


Figure 4.- Comparison of Trailblazer II second-stage spin histories 
computed by tabular , mid -damping, and asymptotic methods with 
more accurate numerical integration results for initial spin 
rate of zero. 5 = 1.75°* 


and to initial spin rate ^i.e., p = + the trial-and-error proc- 

ess is fast and simple after once computing the proportionality factors. 

As might he expected from consideration of the errors involved, the asymp- 
totic method appears to he slightly more accurate than either of the other two 
approximation methods; however, all approximations are within about 2 percent 
of the numerical solution. 

In an effort to exploit the weakness of the mid-damping method, problem 1 
was rerun with an initial spin rate of 200 rad/sec, a value considerably above 
steady-state roll. These results, presented in figure 5* show that all 
approximation methods again compare quite well with the numerical solution. 

As in problem 1, the spin-rate error generated by assuming exponential depend- 
ence of the driving and damping terms is illustrated by the difference between 
the asymptotic and numerical solutions. 


15 



pin rate, rad/sec 



CONCLUDING REMARKS 


An approximate analytical method (with three alternate integration pro- 
cedures) is developed for solving the single -degree -of -freedom roll equation 
with time-dependent coefficients. The method is applied with each integration 
procedure to two sample problems and found to compare closely with more exact 
numerical solutions. The closed-form solution avoids tedious step-by-step 
integration and allows rapid hand computation of results. Of the three inte- 
gration procedures presented, the asymptotic method, if convergent, was 
judged superior for solving the roll equation. Information governing the 
approximate error of each integration procedure is also presented. 


Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., August 13, 1965- 


17 



APPENDIX 


ERRORS OF APPROXIMATION 


Tabular Method 

Errors considered herein are those resulting from the approximation 


eat « i + at + 

2 

The first-order approximation of the general error expression for p is 


where 


Ap = _§E_ Ae at = dp Ae at 
at -3 at 

be 4e 


AfiQ-t — e at _ 



(Al) 


dp 


p _ Be^t + Ce 0 ^ - p(Ae a_ t + d) 


de' 


at 


ae 


at 


ae 


at 


Mid-Damping Method 

The use of a constant (midinterval) damping coefficient over an interval 
results in errors determined approximately as follows. No curve-fitting errors 
were considered. Equation (4) without the C term can be written as 


where 


p = -Aj_p + Be^t* 


A 1 = Ae at + D 


(A2) 
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APPENDIX 


I he first-order approximation of the general error expression for p is 


Ap = AA]_ + -E. Ap + — — — — 
BAi bp S(Be bt ) 


^Be bt ) 


(A3) 


and from equation (A2) 



However, for best results, these slopes should be evaluated at their mean value 
over the associated error increment. Thus, 


Also, 




dt 



"Al,meanAP 




dt 




Since in equation (A3), Ap 
equation can be written as 


and A^Be^) = 0, the approximate error 
dt x ' 


a(A P ) 

dt 



paAe 8 ^ 


dt 



Ae at £ 
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APPENDIX 


where values of p and p are obtained from equation (ll). Finally, by 
integration, 

f t d(Ap) k 2 F (a-Ai)t 1 

Ap = / dt = kit + e' - 1 

J 0 ^ A X - a L J 

- ^rE (a+b)t - d ] - 'f* ( p - ») dt {kk) 

where 

k x = k 2 e( a_Al ) t + kje ( a+b ) b 



=44 

Ai + b) 
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